capture program drop ccn_het_newsym_flexdelta
program define ccn_het_newsym_flexdelta, eclass
	local model_type = "`1'"
	local k = "`2'"
	local inv_var = "`3'"
	local dep_var = "`4'"
	local se_type = "`5'"
	local controls = "`6'"	
	local swap_type = "`7'"
	local beta_het = "`8'"
	local beta_controls = "`9'"
	local delta_controls = "`10'"
	local Iinter = "`11'"
	preserve
	local control_length = $control_length
	capture qui rename `inv_var' I
	capture qui rename `dep_var' S
	qui gen cens_ind = (I==0)
	qui sum I
	local samp = r(N)
	*will need a variable where we drop the I=2080
	*note, I use Imod for I in all cases to compute imr term, across all methods.
	qui gen Imod = I
	qui replace Imod = . if Imod == 2080
	gen cens_indmod = (Imod == 0)
	if ("`model_type'" == "naiveN") {
		qui sum I
		local samp = r(N)
		if (`beta_het' == 0) {
			if ("`se_type'" == "simple") {
				qui reg S I
			}		
			matrix COEFFS = e(b)
			ereturn post COEFFS
			ereturn scalar BETA = _b[I]
			ereturn scalar DELTA = .	
			ereturn scalar DELTA_SE = .
		}
		if (`beta_het' == 1) {
			if ("`se_type'" == "simple") {
				*need to create I'Z
				foreach control_var in `beta_controls' {
					gen Ibeta_`control_var' = I*`control_var'
				}
				reg S I Ibeta_*
			}
			matrix COEFFS = e(b)
			matrix VARS = e(V)
			ereturn post `b' `V'
		}

		di "inputs used: K=`k'; inv_var=`inv_var'; dep_var=`dep_var'; se_type=`se_type', beta_het = `beta_het'"
	}	
	if ("`model_type'" == "naive") {
		qui sum I
		local samp = r(N)
		if (`beta_het' == 0) {
			if ("`se_type'" == "simple") {
				reg S I dumX_`k'_* `controls'
				matrix COEFFS = e(b)
				ereturn post COEFFS
				ereturn scalar ESAMP = e(N)
				ereturn scalar SAMP = `samp'		
				ereturn scalar DELTA = .	
				ereturn scalar BETA = _b[I]
			}
		}
		if (`beta_het' == 1) {
			*need to create I'Z
			di `beta_controls'
			foreach control_var in `beta_controls' {
				gen Ibeta_`control_var' = I*`control_var'
			}
			if ("`se_type'" == "simple") {
				if (`Iinter' == 1) {
					reg S Ibeta_* dumX_`k'_* `controls'
				}
				if (`Iinter' == 0) {
					reg S I Ibeta_* dumX_`k'_* `controls'
				}
				matrix COEFFS = e(b)
				matrix VARS = e(V)
				ereturn post `b' `V'
			}
		}
	}

	if ("`model_type'"=="het_uniform") {
		qui gen imr_term = .
		levelsof X_`k', local(levels)
		foreach i of local levels {
			qui sum Imod if Imod>0 & X_`k' == `i'
			local pos_mean = r(mean)
			qui sum cens_indmod if  X_`k' == `i'
			local bunching = r(mean)
			qui replace imr_term = -`pos_mean'*(`bunching'/(1-`bunching')) if  X_`k' == `i'
		}
		qui gen reg_term = imr_term*cens_ind + I
		foreach i of local levels {
			qui gen het_reg_term_`i' = 0
			qui replace het_reg_term_`i' = reg_term if X_`k' == `i'
		}
		if (`beta_het' == 0) {
			if ("`se_type'" == "simple") {
				reg S I dumX_`k'_* `controls' het_reg_term_*, noconstant
				ereturn scalar ESAMP = e(N)
				ereturn scalar SAMP = `samp'		
				ereturn scalar DELTA = _b[reg_term]	
				ereturn scalar DELTA_SE = _se[reg_term]
				ereturn scalar BETA = _b[I]
				ereturn scalar SE = _se[I]		
			}
		}
		if (`beta_het' == 1) {
			if ("`se_type'" == "simple") {
				*need to create I'Z
				foreach control_var in `beta_controls' {
					gen Ibeta_`control_var' = I*`control_var'
				}
				foreach control_var in `delta_controls' {
					gen reg_term_delta_`control_var' = reg_term*`control_var'
				}
				if (`Iinter' == 1) {
					reg S reg_term reg_term_delta_* Ibeta_* dumX_`k'_* `controls' het_reg_term_* , noconstant
				}
				if (`Iinter' == 0) {
					reg S reg_term reg_term_delta_* I Ibeta_* dumX_`k'_* `controls' het_reg_term_* , noconstant
				}
				matrix COEFFS = e(b)
				matrix VARS = e(V)
				ereturn post `b' `V'	
			}
		}
		di "inputs used: K=`k'; inv_var=`inv_var'; dep_var=`dep_var'; se_type=`se_type', beta_het = `beta_het'"
	}

	if ("`model_type'" == "tobit") {
		qui tobit Imod dumX_`k'_*, ll(0) noconstant
		matrix coeff_mat = e(b)
		matrix var_mat = e(V)
		qui predict imr_term, e(.,0)
		local sigma = coeff_mat[1,`k'+2]^0.5

		qui gen reg_term = imr_term*cens_ind + I
		qui levelsof X_`k', local(levels)
		foreach i of local levels {
			qui gen het_reg_term_`i' = 0
			qui replace het_reg_term_`i' = reg_term if X_`k' == `i'
		}
		foreach i of local levels {
			qui sum imr_term if X_`k' == `i'
			local imr_term_`i' = r(mean)
		}
		if (`beta_het' == 0) {
			if ("`se_type'" == "simple") {
				reg S I dumX_`k'_* `controls' het_reg_term_* , noconstant
				ereturn scalar ESAMP = e(N)
				ereturn scalar SAMP = `samp'		
				*ereturn scalar DELTA = _b[reg_term]	
				*ereturn scalar DELTA_SE = _se[reg_term]
				ereturn scalar BETA = _b[I]
				ereturn scalar SE = _se[I]		
			}
		}
		if (`beta_het' == 1) {
			if ("`se_type'" == "simple") {
				*need to create I'Z
				foreach control_var in `beta_controls' {
					gen Ibeta_`control_var' = I*`control_var'
				}
				foreach control_var in `delta_controls' {
					gen reg_term_delta_`control_var' = reg_term*`control_var'
				}
				if (`Iinter' == 1) {
					reg S reg_term reg_term_delta_* Ibeta_* dumX_`k'_* `controls' het_reg_term_*, noconstant
				}
				if (`Iinter' == 0) {
					reg S reg_term reg_term_delta_* I Ibeta_* dumX_`k'_* `controls' het_reg_term_*, noconstant
				}
				matrix COEFFS = e(b)
				matrix VARS = e(V)
				ereturn post `b' `V'	
			}
		}
		di "inputs used: K=`k'; inv_var=`inv_var'; dep_var=`dep_var'; se_type=`se_type', beta_het = `beta_het'"
	}
	if ("`model_type'" == "het_tobit") {
		*Run Tobit
		qui gen imr_term = .
		qui levelsof X_`k', local(levels)
		foreach i of local levels {
			qui capture tobit Imod if X_`k' == `i', ll(0)
			if _rc==0 {
				qui matrix coeff_mat_`i' = e(b)
				qui matrix var_mat_`i' = e(V)
				qui predict imr_temp, e(.,0)
				qui replace imr_term = imr_temp if X_`k' == `i'
				qui drop imr_temp
				local sigma_`i' = coeff_mat_`i'[1,2]^0.5
			}
			else {
				drop if X_`k'==`i'
			}
		}
		*get an imr_term (cens_exp) for each keep
		qui levelsof X_`k', local(levels)
		foreach i of local levels {
			qui sum imr_term if X_`k' == `i'
			local imr_term_`i' = r(mean)
		}
		qui gen reg_term = imr_term*cens_ind + I
		foreach i of local levels {
			qui gen het_reg_term_`i' = 0
			qui replace het_reg_term_`i' = reg_term if X_`k' == `i'
		}

		********************************************************************************
		*Corrected Models
		if (`beta_het' == 0) {
			if ("`se_type'" == "simple") {
				qui reg S I dumX_`k'_* `controls' het_reg_term_*		
				ereturn scalar ESAMP = e(N)
				ereturn scalar SAMP = `samp'		
				*ereturn scalar DELTA = _b[reg_term]	
				*ereturn scalar DELTA_SE = _se[reg_term]
				ereturn scalar BETA = _b[I]
				ereturn scalar SE = _se[I]		
			}
		}
		if (`beta_het' == 1) {
			*need to create I'Z
			foreach control_var in `beta_controls' {
				gen Ibeta_`control_var' = I*`control_var'
			}
			foreach control_var in `delta_controls' {
					gen reg_term_delta_`control_var' = reg_term*`control_var'
				}
			if ("`se_type'" == "simple") {
				if (`Iinter' == 1) {
					reg S reg_term reg_term_delta_* Ibeta_* dumX_`k'_* `controls' het_reg_term_*
				}
				if (`Iinter' == 0) {
					reg S reg_term reg_term_delta_* I Ibeta_* dumX_`k'_* `controls' het_reg_term_*
				}
				matrix COEFFS = e(b)
				matrix VARS = e(V)
				ereturn post `b' `V'

			}
		}
		di "inputs used: K=`K'; inv_var=`inv_var'; dep_var=`dep_var'; minlevel=`minlevel'; Qnum = `Qnum'; control_type=`control_type'; Effective sample size = `eff_samp'; Sample size = `samp', beta_het=`beta_het'"
	}

	if ("`model_type'" == "symmetric") {
		*generate a variable that tells whether cluster is censored more than 50%
		levelsof X_`k', local(levels)
		foreach i of local levels {
			local switch_`i' = 0
			qui sum cens_indmod if X_`k' == `i', d
			if (r(p50) == 1) {
				local switch_`i' = 1
			}
		}
		*Get the Censored Expectations
		*Step 1 : get the censorted expectations assuming symmetry.
		*These will not be correct within a cluster if the censoring is more than 50%
		gen cens_exp = .
		foreach i of local levels {
			*Step 1 -- find F_{X|Z=z)(0) -- this is just the proportion of observations with Z=z who have X_i = 0
			qui sum cens_indmod if X_`k' == `i'
			*qui gen hatF_0_`i' = r(mean)
			qui local hatF_0_`i' = r(mean)
			qui local op_hatF_0_`i' = 1- `hatF_0_`i''

			*Step 2: find the 1-hatF_0_i percentile amon X (I) such that Z = i
			qui cumul Imod if X_`k' == `i', gen(Icumul_`i') equal
			qui sum Icumul_`i' if X_`k' == `i'
			local Icumul_min = r(min)
			if `op_hatF_0_`i'' > `Icumul_min' { 
				qui sum Imod if X_`k' == `i' & Icumul_`i' <= `op_hatF_0_`i''
				local step2_`i' = r(max)
			}
			if `op_hatF_0_`i'' <= `Icumul_min' {
				local step2_`i' = 0
			}
			*Step 3: find the average value of I conditional on X >= step2 and Z = i
			qui sum Imod if Imod >= `step2_`i'' & X_`k' == `i'
			local step3_`i' = r(mean)

			*Step 4: Calculate the key censored expectation, conditional on Z = i
			local cens_expec_`i' = `step2_`i''- `step3_`i''
			qui replace cens_exp = `cens_expec_`i'' if X_`k' == `i'

		}
		*also want switch codes for cases where the below qreg does not work.
		foreach i of local levels {
			qui capture qreg Imod if X_`k' == `i', q(`op_hatF_0_`i'')
			if (_rc != 0) {
				local switch_`i' = 1
				local qreg_switch_`i' = 1
			}
		}
		*get total number of switches
		local tot_switch = 0
		foreach i of local levels {
			local tot_switch = `tot_switch' + `switch_`i''
		}
		ereturn scalar CLUS_SWITCH = `tot_switch'
		if (`tot_switch'>0) {
			if ("`swap_type'" == "tobit") {
				*Now use standard Tobit
				qui tobit Imod dumX_`k'_*, ll(0) noconstant
				matrix coeff_mat_tob = e(b)
				matrix var_mat_tob = e(V)
				qui predict imr_term_tob, e(.,0)
				local sigma = coeff_mat_tob[1,`k'+2]^0.5

				*get an imr_term (cens_exp) for each keep
				foreach i of local levels {
					qui sum imr_term if X_`k' == `i'
					local imr_term_`i' = r(mean)
					*swap in if the switch term is on
					if (`switch_`i'' == 1) {
						local cens_exp_`i' = `imr_term_`i''
						replace cens_exp = `imr_term_`i'' if X_`k' == `i'
					}
				}

			}
			if ("`swap_type'" == "het_tobit") {
				*Run Tobit
				qui gen imr_term = .
				foreach i of local levels {
					if (`switch_`i'' == 1) {
						qui capture tobit Imod if X_`k' == `i', ll(0)
						if _rc==0 {
							qui matrix coeff_mat_`i' = e(b)
							qui matrix var_mat_`i' = e(V)
							qui predict imr_temp, e(.,0)
							qui replace imr_term = imr_temp if X_`k' == `i'
							qui drop imr_temp
							local sigma_`i' = coeff_mat_`i'[1,2]^0.5
						}
						else {
							drop if X_`k'==`i'
						}
						*get an imr_term (cens_exp) for each keep
						qui sum imr_term if X_`k' == `i'
						local imr_term_`i' = r(mean)
						replace cens_exp = `imr_term_`i'' if X_`k' == `i' //asign tobit piece to relevant bucket.
					}
				}
			}
		}
		qui gen reg_term = cens_exp*cens_ind + I
		levelsof X_`k', local(levels)
		foreach i of local levels {
			qui gen het_reg_term_`i' = 0
			qui replace het_reg_term_`i' = reg_term if X_`k' == `i'
		}
		noisily sum het_reg_term_*
		
		********************************************************************************
		*Corrected Models
		if (`beta_het' == 0) {
			if ("`se_type'" == "simple") {
				qui reg S I dumX_`k'_* `controls' het_reg_term_*		
				ereturn scalar ESAMP = e(N)
				ereturn scalar SAMP = `samp'		
				*ereturn scalar DELTA = _b[reg_term]	
				*ereturn scalar DELTA_SE = _se[reg_term]
				ereturn scalar BETA = _b[I]
				ereturn scalar SE = _se[I]		
			}
		}
		if (`beta_het' == 1) {
			*need to create I'Z
			foreach control_var in `beta_controls' {
				gen Ibeta_`control_var' = I*`control_var'
			}
			foreach control_var in `delta_controls' {
					gen reg_term_delta_`control_var' = reg_term*`control_var'
				}
			if ("`se_type'" == "simple") {
				if (`Iinter' == 1) {
					reg S reg_term reg_term_delta_* Ibeta_* dumX_`k'_* `controls' het_reg_term_*
				}
				if (`Iinter' == 0) {
					reg S reg_term reg_term_delta_* I Ibeta_* dumX_`k'_* `controls' het_reg_term_*
				}
				matrix COEFFS = e(b)
				matrix VARS = e(V)
				ereturn post `b' `V'

			}
		}
	}
	restore
end
